function [likel,ssvec,flag_ok,lht]=...
likelNanSplitsGeneral(parest,parStru,posStru,model,dataStru,filterStru,flags)
% ========================================================================= 
% function [likel,ssvec,flag_ok,lht]=...
% likelNSplitsGeneral(parest,parStru,posStru,model,dataStru,filterStru,flags)
%
%% LikelNSplitsGeneral.m
%
% This code estimates the likelihood of a DSGE model with potentially 
% 1) NaNs in the observables  model 
% 2) NB Breaks, where the solution matrices have NB pages 
% 
% No temporal aggregation 
% 
%% Inputs 
% parest:       portion of parameter vector being estimated
%
% parStru:      structure with
%               .param if the model is solved using a vector of parameters
%           OR  .namesEstimated if the model is solved using
%
% posStru:      .param.est  rows of paramter vector estimated
%
% model:        .handle
%               .addsol
%               .solveOpt
%
% dataStru:     .data     [T n] vector of data
%               .trainVec [2 1] position to start and end the sample in
%                         computing the likelihood
%
% filterStru:   .tauVec   vector of indicators of which pages to use for
%                         each observation
%               .aZero    Initial state,  if flags.initialCond==1
%               .pZero    Initial VCV,    if flags.initialCond==1
%
% flags        .initialCond
% =========================================================================
%% Initialize parameters
likel=-1e20; flag_ok=0;
%% Keep this always to ensure working with latest Value
paramVec=parStru.param; 
paramVec(posStru.param.est)=parest;
%% Model solution, second sample stored in structure second (second sample)
[G,R,C,eu,SDX,Z,structOne,ssvec,structTwo]=feval(model.handle,paramVec,model.solveOptions,model.addsol);
if ~isequal(eu,[1;1])
    lht=[]; 
    return 
end 

[T,ny]=size(dataStru.data);
dataStru.data = dataStru.data';
%% Initialization
%for j = 1:15
if flags.initialCond==0
    MM=R(:,:,filterStru.tauVec(1))*(SDX(:,:,filterStru.tauVec(1))');
    pshat=feval(@lyapunov_symm,G(:,:,1),MM*(MM'));
    %     pshat=feval(@lyapunov_symm,G(:,:,filterStru.tauVec(1)),...
    %         R(:,:,filterStru.tauVec(1))*(SDX(:,:,filterStru.tauVec(1))')*...
    %         SDX(:,:,filterStru.tauVec(1))*(R(:,:,filterStru.tauVec(1))'));
    shat       =zeros(size(G,1),1);
else
    pshat=filterStru.pZero;
    shat=filterStru.aZero;
end
lht=zeros(T,1);
W          = eye(ny);
Zdim       = zeros(T,1);
% Filter

for ii=1:T;
    % Handling of missing observations 
    ytt=dataStru.data(:,ii);    
    % Determine W and position of the NAN  
    ind =~isnan(ytt);    
    ytt=ytt(ind);    
    Zdim(ii)=length(ytt); 
    Ztt=W((ind==1),:)*Z(:,:,filterStru.tauVec(ii));
    ytt=ytt-Ztt*C(:,filterStru.tauVec(ii));    
    [shat,pshat,lht(ii)]=feval(@kf,ytt,...
        Ztt,shat,pshat,G(:,:,filterStru.tauVec(ii)),R(:,:,filterStru.tauVec(ii))*(SDX(:,:,filterStru.tauVec(ii))'));
end

%{
end

tempVEC = dataStru.start:0.25:dataStru.finish;
tempVECind = tempVEC>=dataStru.startest;

for j = 1:15
    subplot(3,5,j);
    plot(yhat{j}(tempVECind))
end
%pause
%}
%% Likelihood plus integration constant
%disp(sprintf('Likel Before Estimation %10.4f',lht(dataStru.trainVec(1)-1))); 
lht   = lht(dataStru.trainVec(1):dataStru.trainVec(2));
Zdim  = Zdim(dataStru.trainVec(1):dataStru.trainVec(2)); 
% Cons1=-0.5*length(lht)*ny*log(2*pi); 
% Cons2= -0.5*(sum(Zdim))*log(2*pi); 
% fprintf('Differences in constants %2.10f \n',Cons1-Cons2); 
likel = -0.5*(sum(Zdim))*log(2*pi)+sum(lht);
flag_ok=1;
%% End of File
end